4

I started out by posting https://stackoverflow.com/questions/6963991/how-can-i-transform-the-coordinates-of-a-shapefile

In short, I found a shapefile here, and I want to be able to convert that to Latitude/Longitude (Geodetic?) coordinates that I can store in my PostGIS database. I am having trouble working with the Ruby libraries, but that is a separate issue. First, I need to correctly understand the projection.

In the linked zip file, there is a .prj file, which contains the following (It is only one line, but I broke it up for readability):

PROJCS["NAD_1983_StatePlane_California_III_FIPS_0403_Feet",
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",
SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",6561666.666666666],PARAMETER["False_Northing",1640416.666666667],
PARAMETER["Central_Meridian",-120.5],PARAMETER["Standard_Parallel_1",37.06666666666667],
PARAMETER["Standard_Parallel_2",38.43333333333333],PARAMETER["Latitude_Of_Origin",36.5],
UNIT["Foot_US",0.3048006096012192]]

I was unable to find tools that translate .prj files into a projection string, so I used this as a reference and came up with:

+proj=lcc +datum=NAD83 +ellps=GRS80 +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=6561666.666666666 +y_0=1640416.666666667 +units=us-ft +pm=greenwich +no_defs

I happen to know that the points in this shapefile should correspond to locations in San Francisco, so I could test my work:

cs2cs -f "%.8f" +proj=lcc +datum=NAD83 +ellps=GRS80 +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=6561666.666666666 +y_0=1640416.666666667 +units=us-ft +pm=greenwich +no_defs +to +proj=lonlat +datum=WGS84 +ellps=WGS84
6011287.4999795845 2100857.2499904726
-164.18117482   17.72087079 0.00000000

That's about 1000km west of Hawaii, so clearly that's not right. I tried many different things, but the closest I got was by stripping the +units=us-ft off.

cs2cs -f "%.8f" +proj=lcc +datum=NAD83 +ellps=GRS80 +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=6561666.666666666 +y_0=1640416.666666667 +pm=greenwich +no_defs +to +proj=lonlat +datum=WGS84 +ellps=WGS84
6011287.4999795845 2100857.2499904726
-126.98864867   40.47509388 0.00000000

That's about 500 km NW of San Francisco.

What am I missing here?
If I import the shapefile into the online ArcGIS Explorer, it understands it perfectly.

Improve this question

3 Answers 3

Reset to default
4

I believe that the correct string for this projection is: 

+proj=lcc +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000 +y_0=500000.0000000002 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs

See: http://spatialreference.org/ref/esri/102643/

Improve this answer
2
  • 1
    Just to clarify, the proj4 addendum ftp.remotesensing.org/proj/proj.4.3.pdf (p. 7) explains that the + parameters always need to be specified in metres, regardless of the units the projection is in. That's why the +x and +y definitions in Nick's answer are different, because Larry's parameters were in us-ft. This is something that has caught me out in the past.
    – MerseyViking
    Commented Aug 10, 2011 at 9:28
  • Yes! This seemed to work. Thank you SO much. And thanks for the explanation MerseyViking.
    – Larry
    Commented Aug 10, 2011 at 15:55
2

Coordinate reference systems and transformations are a complicated matter and are a source of trouble in practice.

I'm not aware of any tool that converts WKT to proj4. As a general hint, look up the PROJCS keyworks on the Spatial Reference website (that's how nhopton found his answer) to find hints to your projection.

Improve this answer
0
2

As requested by EPSG31468, here's a very simple Python script to convert from any projection format handled by GDAL to Proj4 strings:

#!/usr/bin/python
# This code is released into the public domain, but may be subject to the licence(s) of the third-party libraries that have been used.
# There is absolutely no warranty or guarantee of fitness for purpose.

import sys
import osgeo.osr as osr

projectionString = sys.stdin.read()

sr = osr.SpatialReference()
sr.SetFromUserInput(projectionString)
proj4 = sr.ExportToProj4()

if proj4 == "":
    print("Unable to convert input string.")
else:
    print(proj4)

There's not a great deal of error handling, but it works with the strings I've tried it with.

To use it, type something like:

$ cat larrysWKT.txt | toproj4.py
+proj=lcc +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000 +y_0=500000.0000000001 +ellps=GRS80 +units=us-ft +no_defs

or even:

$ echo EPSG:27700 | toproj4.py
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs

Equivalently, Windows users will need to do this:

c:\> type larrysWKT.txt | python toproj4.py
c:\> echo EPSG:27700 | python toproj4.py

Naturally you'll need Python installed, as well as GDAL/OGR and the corresponding Python bindings. Linux users should be able to get this from their package repositories, and I think Windows users can install it from the OSGeo4W installer.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.